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Introduction 


In  situ  measurement  of  thin  film  growth  is  of  great  interest  to  both  scientists  and 
engineers.  For  the  scientist,  real-time  measurement  of  film  properties  within  the 
deposition  chamber  enables  the  study  of  evolution  and  growth  processes.  For  the 
engineer,  such  techniques  offer  greater  process  control  and,  ultimately,  improve  devices. 
Among  those  methods  which  can  adapted  for  in  situ  measurement,  x-ray  reflectivity  is 
particularly  well  suited  to  the  characterization  of  thin  films.  X-ray  reflectivity  can  recover 
density,  thickness,  and  roughness  information  from  layers  and  multilayers  through  bulk 
reflection  and  interference  effects.  It  is  equally  capable  of  measuring  metals  and 
dielectrics,  unlike  ellipsometry,  and  the  short  wavelength  of  the  probing  radiation  enables 
the  measurement  of  films  having  thicknesses  as  low  as  a  few  nanometers  [1].  In  addition, 
x-rays  are  insensitive  to  environment,  making  x-ray  reflectivity  suitable  for  use  in 
sputtering  chambers,  unlike  electron-beam  techniques  [2], 

Conventionally,  x-ray  reflectivity  measurements  have  probed  k-space  by  varying 
the  incident  angle  of  a  monochromatic  beam  in  an  angular  dispersive  technique.  This 
methodology  has  been  difficult  to  integrate  into  deposition  systems  by  virtue  of  the 
expense  and  time  delay  inherent  to  high  precision  goniometer  scans.  The  system 
developed  as  part  of  this  study  circumvents  these  concerns  by  utilizing  an  energy 
dispersive  technique.  A  fixed  beam  geometry  is  employed  and  k-space  probed  by 
photons  having  a  range  of  incident  energies.  Recent  work  by  Chason  [3]  and  Kellerman 
[4]  has  shown  the  suitability  of  this  technique  for  use  in  situ.  At  present,  work  done  on 
the  measurement  of  growing  films  has  concentrated  on  the  use  of  single  energies  and 
fixed  angles  to  produce  “deposition  curves”-  variations  in  reflectivity  with  film 
growth[5,6,7].  Energy  dispersive  technique  expands  on  these  methods  to  produce  a 
standard  reflectivity  curve  which  evolves  in  time. 

In  order  to  study  the  use  of  energy  dispersive  x-ray  reflectivity  for  in  situ 
characterization,  a  specialized  deposition  and  measurement  system  has  been  developed. 
The  system  makes  use  of  an  existing  x-ray  diffractometer  through  a  demountable  vacuum 
deposition  chamber,  computer  based  acquisition  card,  and  custom  analysis  software. 
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Theory  of  X-ray  Reflectivity 


Optical  means  of  material  characterization  differ  from  one  another  in  two  aspects: 
the  first  is  the  energy  of  the  incident  photons  used  and  the  second  is  the  geometry  chosen 
for  the  positioning  of  source  and  detector.  In  the  technique  of  x-ray  reflectivity, 
information  is  gathered  by  photons  having  energies  in  the  range  of  1  to  100  KeV  that  are 
directed  at  the  material  at  small  grazing  angle  and  detected  in  the  region  of  specular 
reflection.  The  interaction  of  materials  with  radiation  in  the  x-ray  region  of  the  spectrum 
is  quite  different  from  that  in  the  visual  region.  Due  to  their  high  frequency,  x-rays  do 
not  interact  with  nuclei  but  instead  are  sensitive  to  electron  density  [8].  Moreover,  the 
high  energies  involved  allow  for  absorption  due  to  electronic  ionization.  Both  of  these 
properties  are  taken  into  account  by  writing  the  index  of  a  material  as  n  =  1  -  6  -  ip.  In 
this  expression,  5  determines  the  phase  change  of  the  wave  and  is  directly  related  to 
electron  density,  while  P  determines  the  x-ray  absorption.  The  two  constants  are  related 
to  the  material  properties  by 

<5  =  ^an<iA  =  A^W 
271  Ati 

where  p  is  the  electron  density,  ro  is  the  classical  electron  radius,  and  |i(X)  is  the  linear 
absorption  coefficient 

Both  5  and  p  are  small  positive  values.  Thus,  for  most  materials,  the  index  of 
refraction  for  x-rays  will  be  slightly  less  than  one.  This  would  appear  problematic  due  to 
the  suggestion  of  a  wave  speed  greater  than  the  speed  of  light  in  vacuum,  c.  However, 
while  the  index  determines  the  phase  velocity  of  the  wave,  the  group  velocity  (the  speed 
at  which  energy  is  transported)  remains  less  than  c.  [8]  Because  the  refractive  index  is 
less  than  one,  x-rays  passing  from  air  into  a  material  can  experience  total  external 
reflection  over  a  range  of  incident  grazing  angles  between  0  and  some  critical  angle,  6c. 
The  critical  angle  can  be  calculated  from  Snell’s  law,  which  is  written  as 

Wj  cos^,  =  cos  ^2 

for  angles  measured  from  the  material  surface.  In  general,  P  is  smaller  than  6  by  an  order 
of  magnitude.  Therefore,  in  calculating  the  critical  angle,  the  index  of  a  material  can  be 
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taken  as  n  =  1  -  6.  Using  the  Taylor  expansion  for  the  cosine,  the  critical  angle  can  be 
related  to  the  electron  density  of  the  material. 


cos^^  =  l-S 

1 — ‘-  =  \-5 
2 

e,  =  ^[26  oc  ^ 


From  this  expression  it  is  clear  that  denser  materials  will  have  a  higher  critical  angle.  For 
incident  angles  less  than  0c  the  incident  and  reflected  intensities  will  be  equal.  At  angles 
greater  than  0c  some  portion  of  the  light  will  enter  into  the  material.  This  relation 
between  the  incident  and  reflected  intensity  is  known  as  the  reflectivity 


R  = 


‘  reflected 


‘  incident 


In  the  simple  case  described  here,  the  reflectivity  plot  for  a  material  would  be  a  step 
function:  equal  to  one  at  angles  below  0c  and  dropping  rapidly  at  angles  larger  than  0c.  In 
actual  materials,  absorption  smoothes  this  transition  causing  a  gradual  drop  in  reflected 
intensity  in  the  region  of  the  critical  angle.  Beyond  this  region,  the  reflectivity  can  be 
expressed  in  terms  of  the  incident  and  critical  angles  by 

f  sin  0 


To  express  the  reflectivity  of  a  material  in  greater  detail,  Fresnel  theory  must  be 
used.  Under  the  assumption  that  the  atomic  granularity  of  materials  relative  to  x-ray 
wavelengths  can  be  neglected  due  to  the  relatively  large  penetration  depth  [9],  Fresnel 
theory  for  homogenous  media  can  be  used  to  relate  the  electric  field  amplitudes  of  the 
incoming  and  reflected  waves.  [10]  The  scenario  is  depicted  \n  Figure  1.  For  the  s 
polarization  case  the  Fresnel  reflection  coefficient  is  written  as 


Ei^  _  Hq  sin  9^  - sin 0^  _  -k^ 

Eq  sin  9q  +  n^  sin  0^  k^  +  k^ 


where  A:  ^  is  the  perpendicular  component  of  the  wave  vector  as  defined  by 


kf  =  .  sin  0  =  -^sin  0. . 
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The  reflectivity  is  calculated  from  the  Fresnel  coefficient  by  taking  the  product  of  the 
coefficient  with  its  complex  conjugate 

R  =  rr\ 

Formal  calculation  in  the  case  of  unpolarized  light  requires  that  the  reflectivity  be 
evaluated  as  the  average  of  the  reflectivities  of  s  and  p  polarization.  However,  for 
calculations  involving  short  wavelengths  and  small  angles,  the  effects  of  polarization  are 
negligible  and  the  s  polarization  equations  will  suffice  for  all  cases.  [11] 

For  the  condition  in  which  x-rays  are  incident  from  air  or  vacuum,  the  reflectivity 
from  a  single  layer  can  be  rewritten  in  terms  of  the  incident  angle  and  the  real  and 
complex  parts  of  the  index  [12].  The  first  step  in  the  this  reduction  is  the  application  of  a 
trigonometric  identity 

«,  sin  0i  =  «,  -y/l-cos^  . 

This  becomes 

=  yjni  -  «o  cos^  $0 

by  application  of  Snell’s  law.  Expanding  the  cosine  and  neglecting  terms  second  order  or 
higher  in  6  and  P  yields 

=  40l-25,-2ip,  . 

Using  this  expression  and  replacing  sinGo  with  Go  by  the  small  angle  approximation,  the 
reflection  coefficient  is  now  written  as 

_El  _e,-^el-25,-2ip, 
d,+4el-25,-2ip, 

From  this  form  of  the  reflection  coefficient,  the  reflectivity  can  be  plotted  as  a  function  of 
the  ratio  P  /  6,  as  shown  in  Figure  2. 

Multilayer  Reflectivity 

In  films  which  are  sufficiently  thin,  interference  fringes  can  be  produced  in  the 
reflectivity  curve  by  the  interaction  of  the  reflections  from  the  air/film  interface  and  the 
film/substrate  interface,  as  first  observed  by  Kiessig  [13]. 


5 


Figure  2.  A  plot  of  reflectivity  versus  the  normalized  glancing  angle, 
0/9c,  as  a  function  of  p/5.  As  the  ratio  increases,  the  drop  in  reflectivity 
in  the  region  of  critical  angle  becomes  progressively  more  gradual. 
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An  example  of  this  is  shown  in  Figure  3.  The  spacing  of  the  fringes  is  directly  related  to 
the  thickness  of  the  film  layer.  This  is  expressed  in  terms  of  the  Bragg  condition. 


which  yields, 


d  = 


TV 


For  the  simple  three  layer  interface  of  air,  film,  and  substrate  depicted  in  Figure 
4,  one  can  define  two  sets  of  Fresnel  coefficients.  One  for  the  air/film  interface  and  a 
second  for  the  film/substrate  interface.  Considering  s  polarization  only,  one  has  for  the 
first  layer 

^ ^  _  ”o  sin  6q  -«i  sin 0,  _  -k^ 

sin  9q  +  «,  sin  6^  +  k^  ’ 


^0,1  ~ 


and  for  the  second  layer 


2rto  sin  9q 


2ki 


n.2  ~ 


«o  sin  9^  +  «,  sin  9^  k^  +  k^  ’ 


^0,1  “ 


k2  -K 


k^  +k,^  ’ 


^12  ~ 


w,  sin  9.  k^ 


n^  sin  9^  k^ 

The  reflected  field  for  the  entire  film  will  be  a  superposition  of  the  reflections  from  both 
boundary  layers  with  the  field  from  the  film/surface  interface  being  phase  shifted  by  an 
amount 

^  =  — «,c/sin0,  =k^d . 

^0 

The  electric  field  amplitudes  are  then  written  as 

and 
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By  rearranging  these  expressions  the  total  reflectivity  of  the  layers  can  be 
expressed  in  terms  of  a  transfer  matrix,  of  the  type  commonly  used  for  multilayer  optical 
coatings.  [14] 


1 _ 

=  M,  1, 

i 

_ 1 

The  matrix  of  transmission  and  reflection  coefficients  may  be  expressed  most  generally 
in  terms  of  the  perpendicular  component  of  the  incident  wave  vector,  .  The  transfer 
matrix  is  Avritten  as  the  product  of  the  Fresnel  coefficient  matrix,  F,  and  the  propagation 
matrix,  P. 

M  =  FP, 


where 


and 


0 

exp(-  ikjd) 


in  each  of  the  layers  can  be  calculated  from  the  incident  wave  vector  and  the  index  in 
that  layer  according  to 


kf  =  sin  0j  =  ^ ^0^ -2S^ -2j/3 .  . 

Aq  Aq 

Films  of  graded  density  can  be  modeled  as  a  series  of  layers  of  varying  index,  as  shown 
in  Figure  5.  Multilayer  reflectivity  is  calculated  using  the  product  of  the  transfer 
matrices  for  each  of  the  of  the  layers 


lKa 


7^+1 

j=i 


A  solution  is  found  by  setting  to  zero  under  the  assumption  that  the  thickness  of  the 
substrate  will  prevent  a  reflection  from  the  lower-most  interface.  It  is  then  convenient  to 
define  reflection  and  transmission  coefficients  for  the  entire  system 
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Interface  1 
Layer  1 


Interface  j 


Layer  N 

Interface  N-f  1 
Substrate 


Figure  5.  A  figure  depicting  a  film  composed  of  N  layers  on  a  substrate. 
Air  is  the  0  layer  and  the  substrate  is  the  N+1  layer. 
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-A/’+l 


/-  =  ^  and  t  =  ^ 

/ro  /ro 

Using  these  the  matrix  expression  can  be  rewritten  as 


=r 


where 


N+\ 

y=i 

The  reflection  and  transmission  for  the  layer  is  written  in  terms  of  the  elements  of  the 
transfer  matrix  for  the  system, 

r=^,R  =  rr\ 

-r-\  *  * 


and 


f  =  — ,  T  =  tt\ 


Reflectivity  from  Imperfect  Interfaces 

In  the  preceding  analysis,  ideal  interfaces  between  the  various  media  have  been 
assumed.  However,  in  actual  materials,  the  transition  between  two  regions  having 
different  densities  will  occur  gradually,  as  shown  in  Figure  6.  The  variation  in  index 
across  an  interface  may  be  expressed  in  terms  of  the  error  function 

n{z)=n^  +{n^  -n,)Erf{z,a),  [15] 
where  the  error  function  is  defined  by 

Erf{z, a)  =  (crV^ )  ' exp(- . 

a  in  the  above  expressions  represents  the  standard  deviation  of  the  Gaussian  distribution 
and  can  be  viewed  as  the  rms  interface  width  of  the  transition  between  ni  and  n2,  as 
Figure  6  illustrates. 

As  a  result  of  these  imperfect  interfaces,  scattering  of  the  incident  and  reflected  x- 
rays  will  occur.  The  result  of  this  scattering  at  buried  interface  layers  is  a  reduction  in  the 
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amplitude  of  the  interference  fringes.  Scattering  from  the  surface  of  the  film  increases 
the  rate  of  fall  off  of  the  reflectivity  curve,  in  addition  to  damping  the  interference 
fringes.  This  is  shown  in  Figure  7.  [16] 

Two  models  for  the  scattering  at  interfaces  can  be  considered.  For  both  cases  the 
magnitude  of  the  reflected  component  is  reduced  due  to  the  imperfect  interface.  In  the 
first  model  this  is  assumed  to  produce  a  global  energy  loss  due  to  scattering.  The  second 
model  assumes  that  no  energy  is  lost  and  that,  rather,  the  transmission  through  each  given 
interface  is  correspondingly  increased  [15]. 

The  loss  of  reflectivity  is  incorporated  in  the  model  through  a  modification  of  the 
Fresnel  reflection  coefficients.  Considered  in  terms  of  scattering  formalism  [17],  the 
correction  is  given  by  the  Fourier  transform  of  the  density  gradient  across  the  interface. 

/?'(z)=  (ctV^)  '  exp(-z^/2cr^) 

p'ifi)  =  J  p'{z)e'^^dz  =  exp(-  4cr^  jo) 


which  is  equivalent  to  the  well  known  expression  for  the  Debye-Waller  factor 

Z),  y  =  (- 1 6;7-  V  ^  sin  ^  /Af ) . 


Due  to  assumptions  implicit  in  scattering  formalism,  this  expression  is  only  valid  for  0 
»  0c.  Nevot  et.  al.  [17]  have  developed  more  general  roughness  correction  factors  for 
interfaces  having  predominately  high  spatial  frequencies, 

A./  =exp(-2cr'A,:^A:|), 

and  for  those  having  predominately  low  spatial  frequencies, 

A.>  =exp(-2o-'(Z:,:"y  . 

For  this  study,  the  factor  for  high  spatial  frequencies  is  used.  It  is  included  in  the 
multilayer  model  by  modifying  the  form  of  the  Fresnel  portion  of  the  transfer  matrix  such 
that  the  roughness  factor  multiplies  the  off-diagonal  terms. 
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Figure  7.  A  plot  of  reflectivity  for  a  series  of  similar  40  nm  films  showing  the 
effect  of  surface  and  interface  roughness.  Curve  A  has  no  surface  or  interface 
roughness.  Curve  B  has  an  interface  width  on  the  surface  of  1  nm.  Curve  C 
has  an  interface  width  of  1  nm  at  the  film/substrate  interface. 
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Assuming  the  roughness  model  which  allows  for  global  energy  loss,  this  is  the  only 
correction  necessary. 


Energy  Dispersive  Reflectivity 


In  the  previous  analysis,  the  reflectivity  of  a  multilayer  has  been  expressed  as  a 
function  of  the  component  of  the  propagation  vector  perpendicular  to  the  film  surface, 

.  As  is  evident  fi’om  the  expression  defined  earlier, 

^  W;  sin  ^ ^  , 


k^  is  a  function  of  both  the  wavelength  of  the  incident  beam,  Xo,  and  the  angle  of 
incidence,  9o.  Conventional  practice  has  followed  after  the  initial  work  of  Kiessig  [13] 
and  Parratt  [9],  both  of  whom  measured  reflectivity  by  scanning  through  0  with  a 
monochromatic  beam  and  plotting  the  intensity  of  the  specular  reflection  measured  at  26. 
However,  it  is  possible  to  gather  the  same  information  by  measuring  the  intensity  as  a 
function  of  wavelength  using  a  fixed  angle,  specular  geometry.  In  this  method,  first  used 
by  Bilderback,  et  al.  [18]  to  analyze  the  critical  angle  of  x-ray  mirrors,  the  incident  and 
outgoing  angles  are  fixed  and  the  energy  spectrum  of  a  white  source  is  profiled.  Dhez,  et. 
al  [19]  have  illustrated  the  relationship  of  these  two  techniques  graphically  \n  Figure  8. 

A  comparison  of  two  reflectivity  curves  of  the  same  sample  made  by  the  angular 
dispersive  and  energy  dispersive  techniques  is  shown  in  Figure  9  [1] . 

While  the  energy  dispersive  technique  allows  for  a  simplified  mechanical 
apparatus,  analysis  of  energy  dispersive  data  presents  additional  challenges.  The  index  of 
materials  for  x-rays  has  been  expressed  in  terms  of  real  and  imaginary  components,  5  and 
P.  Each  of  these  has  a  dependence  on  the  wavelength  of  the  x-rays.  The  behavior  of  5 
and  p  can  be  predicted  by  the  classical  dispersion  laws 
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Figure  8.  Graph  illustrating  the  relationship  between  energy  dispersive  and 
angular  dispersive  methods  of  measurement.  The  reflectivity  of  a  periodic 
multilayer  is  depicted  as  a  function  of  both  angle  and  energy.  The  maximum 
intensity  of  a  given  order  of  Bragg  peak  will  always  lie  along  the  line  drawn 
in  the  horizontal  plane. 
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Figure  9.  Angular-dispersive  reflectivity  compared  with  fixed  angle 
results  for  the  three  tantalum  on  silicon  samples.  Curve  (a)  gives  normal 
reflectivity  data  (solid  curve)  compared  with  0.6°  20  fixed  angle  data 
(closed  triangles)  and  1.5°  20  fixed  angle  data  (open  triangles)  for  sample 
1 .  Curve  (b)  gives  the  normal  reflectivity  data  (solid  curve)  compared 
with  0.6°  20  fixed  angle  (closed  squares)  and  1.5°  20  fixed  angles  (open 
squares)  for  sample  2.  Curve  (c)  shows  the  normal  reflectivity  (solid 
curve)  versus  the  0.6°  20  fixed  angle  (closed  circles)  and  the  1.5°  20  fixed 
angle  (open  circles)  information  for  sample  3.  By  varying  the  incident 
angle,  the  limited  range  of  incident  energies  can  be  used  to  profile  the 
entire  reflectivity  curve.  Often,  however,  only  a  range  of  the  curve  is 
necessary  for  analysis. 
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where  f  is  the  real  dispersion  correction  (due  to  the  scattering  cross-section)  and  f’  is 
the  imaginary  dispersion  correction.  However,  extrapolating  the  index  based  on  these 
formulas  is  only  effective  for  small  changes  in  X.  The  indices  may  also  may  computed 
from  f  and  f  ’  listed  in  tables;  however,  current  tables  do  not  extend  beyond  10  keV.  In 
practice  the  values  are  predicted  from  theory,  but  fit  experimentally  as  two  additional  free 
parameters. 

Reflectivity  from  Growing  Films 


For  a  growing  film  thickness,  density,  and  surface  roughness  will  be  evolving  in 
time.  Modeling  such  a  system  as  a  single  layer  bound  by  atmosphere  and  substrate 
layers,  as  shown  in  Figure  4,  a  total  reflection  coefficient  is  written  as, 

Ej  1  +  ''o.i''i.2^o.iA.2  exp(2r<ife/- ) 

For  regions  of  the  curve  far  from  the  critical  angle  the  Fresnel  coefficients  take  on  the 
asymptotic  form 


w. 


Asin^e 


[5]. 


If  the  average  index  of  the  film  is  assumed  to  be  constant,  the  reflectivity  of  the  film  can 
be  considered  a  function  of  d,  a,  and  .  Reflectivity  measured  using  fixed  angle 
geometry  will  evolve  along  with  d  and  a.  As  the  layer  thickens,  the  reflectivity  will  vary 
sinusoidally  as  the  wave  reflected  from  the  film  surface  and  the  wave  reflected  from  the 
substrate  pass  in  and  out  of  the  constructive  interference  condition.  Assuming  that  3  can 
be  neglected  in  the  expression  for  the  index,  the  reflectivity  can  be  expressed  in  terms  of 
the  Bragg  relation 


m 


R{t)=2T  1  — 

L  si 


sin  ^ 


[5], 


where  T  is  the  spatial  period  of  the  oscillations  and  m  is  order  of  interference.  As  the 
film  thickens  it  roughens,  damping  the  amplitude  of  the  peaks  in  the  reflectivity  plot. 
This  growth  in  roughness  can  be  modeled  as 
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w{d)  =  t^  [201 

where  3  ranges  between  .2  and  1.  Baranov  et.  al,  [5],  Lee  et.  al.  [6],  and  Heilman  et.  al. 
[7]  have  demonstrated  this  technique  for  evaporated  and  sputtered  films.  An 
experimental  curve  for  the  growth  of  aluminum  by  evaporation  is  shown  in  Figure  10. 

A  primary  concern  of  using  x-ray  reflectivity  for  the  measurement  of  growing  thin 
film  is  the  time  averaging  which  will  occur  as  a  result  of  the  finite  collection  time 
required  to  generate  a  cure.  The  reflectivity  of  a  sample  will  vary  as  it  thickens  and 
roughens.  For  a  collected  curve  to  be  meaningful,  the  variation  in  the  reflectivity,  AR, 
must  be  small  over  the  collection  time,  At,  relative  to  the  features  of  the  curve  which  are 
of  interest.  For  this  reason,  angular  dispersive  technique  is  not  feasible.  The  fixed  angle 
technique  described  is  far  more  feasible  because  only  x-ray  intensity  and  detector 
efficiency  extend  the  measurement  time.  Nonetheless,  it  is  desirable  to  produce 
conventional  reflectivity  curves  of  the  growing  film.  Using  energy  dispersive  methods,  a 
reflectivity  curve  could  be  constructed  for  each  data  point  on  the  curve  of  a  fixed  angle 
growth  oscillation.  However,  what  must  be  shown  is  that  AR  will  be  sufficiently  small 
over  At  for  useful  curves  to  be  generated.  In  recent  experiments,  a  measurement  time  of 
ten  seconds  has  been  shown  to  produce  a  usable  curve,  though  300  seconds  is  preferred 
for  maximum  resolution  as  shown  in  Figure  11.  [1].  To  investigate  this  a  series  of 
models  were  created  using  spreadsheet  based  x-ray  reflectivity  simulation  software.  A 
layer  of  tantalum  on  a  silicon  wafer  was  considered  with  an  interface  width  between  film 
and  substrate  of  5  nm.  The  average  rms  percentage  errors  between  the  log  reflectivities 
of  the  initial  and  final  states  of  various  growth  conditions  were  calculated.  A  typical 
growth  rate  for  work  with  growth  oscillations  is  .2  nm  min’*  [6] .  For  a  growth  of  1.5  nm 
(450  seconds  measurement  time  at  this  growth  rate)  and  an  increase  in  roughness  given 
by  3  ~  1  of  a  4.0  nm  film  with  initial  interface  width  of  3.0  nm  the  average  error  was  .2%. 
For  a  growth  of  .5  nm  (150  seconds)  the  error  reduced  to  .07%.  Because  thickness 
fiinges  are  often  a  very  small  effect,  a  slow  growth  time  coupled  with  a  fast  measurement 
time  may  be  n^essary.  A  .  1  nanometer  growth  during  a  thirty  second  measurement  time 
results  in  an  error  of  only  .01%  which  is  largely  indistinquishable  to  the  eye  in  the 
semilog  plot  shown  m  Figure  12.  Using  the  fastest  collection  time  often  seconds,  this 
allows  for  minimal  time  smearing  effects  at  growth  rates  up  to  .6  nm  min’*.  In  samples 
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Figure  10.  Growth  oscillation  for  evaporated  aluminum.  The  circles  indicate 
the  reflectivity  measured  as  a  function  of  the  film  thickness.  The  solid  line 
denotes  reflectivity  calculated  from  a  three-layer  model  that  is  assumed  to  have 
surface  roughness  proportional  to  the  film  thickness. 
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Figure  11.  Fixed-angle  energy  dispersive  x-ray  reflectivity  for  a  thin 
tantalum  film  (-8.5  nm)  on  silicon.  The  scan  was  taken  with  20  =  1.5° 
with  copper  source  at  20  kV  and  20  mA.  Curve  (a)  represents  the  energy 
spectrum  of  the  primary  beam  at  29  =  0°  for  300  s.  Curve  (b)  is  the  raw 
energy  data  taken  at  20  =  1.5°  for  300  s.  Curve  (c)  is  the  normalized  data, 
(b)  divided  by  (a),  showing  the  "Kiessig  curves"  or  thickness  oscillations 
for  the  tantalum  film:  triangles  showing  300  s  of  collection  time,  diamonds 
10  s  of  collection  time. 
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Figure  12.  A  theoretical  plot  of  the  reflectivity  curves  of  two  samples 
differing  in  thickness  and  roughness.  This  is  representative  of  a  growth 
time  of  30  seconds  for  a  growth  rate  of  .2  nm  min'\  The  average  rms 
error  between  the  two  curves  is  .01%. 
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having  less  of  a  difference  between  the  indices  of  the  film  and  substrate  the  amplitude  of 
the  thickness  oscillations  will  be  reduced  requiring  lower  deposition  rates. 

Data  Analysis 

Information  may  be  derived  from  a  reflectivity  curve  by  a  number  of  methods. 

As  has  been  described,  the  particular  features  of  a  curve  are  related  to  particular 
properties  of  the  film  which  produced  that  curve.  One  class  of  methods  uses  this  fact  to 
rapidly  determine  a  single  film  property  by  evaluating  a  single  aspect  of  the  reflectivity 
curve.  A  second  class  of  methods  evaluates  all  of  the  film  properties  by  minimizing  the 
fitting  error  between  a  complete  model  of  the  reflectivity  curve  and  the  data  from  the 
experimental  curve.  Both  forms  of  analysis  will  be  used  in  this  study.  For  rapid 
assessment  of  film  thickness  a  Fourier  transform  method  is  used.  More  complete  analysis 
is  made  through  use  of  a  least-squares  fit  to  the  transfer  matrix  model. 

Because  layers  of  different  density  can  give  rise  to  oscillations  in  the  reflectivity 
curve  whose  period  is  related  to  the  thickness  of  those  layers,  it  can  be  expected  that  the 
Fourier  transform  can  be  used  to  recover  thickness  information.  However,  the  intensity 
variations  and  imperfect  periodicity  inherent  in  a  reflectivity  curve  complicate  this 
analysis.  Therefore,  data  must  undergo  a  preliminary  transformation  to  restore  its 
periodicity. 

For  purposes  of  illustration,  the  reflectivity  curve  may  be  represented  as 

where  A  and  B  are  components  which  vary  slowly  with  .  For  a  transform  of  the  data 
to  yield  accurate  results  the  nonperiodic  components  of  the  reflectivity  must  be  removed. 
The  plateau  before  the  critical  angle  and  constant  slope  following,  represented  by  A,  can 
be  removed  by  limiting  the  region  analyzed  and  dividing  by  an  appropriate  factor, 
respectively.  The  varying  periodicity  which  results  fi-om  the  inclusion  of  in  the 

cosine  can  also  be  removed  by  an  appropriate  choice  of  multiplying  factor. 

Scattering  theory  confirms  the  relation  between  the  Fourier  transform  and 
thickness  information.  Under  this  formalism,  the  reflectivity  above  the  critical  angle  of  a 
layer  having  index  profile  n(z)  is  given  in  terms  of  the  scattering  vector,  Q,  by 
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where  the  scattering  vector  represents  the  momentum  transfer  to  the  surface  and  is  related 
to  by 


Q  =  Kn-Ku,='^k^ 


Arc  sin  $ 

I 


From  this  form  of  the  reflectivity  it  is  clear  that  in  the  region  beyond  the  critical  angle  the 
reflectivity  is  proportional  to(A:'‘'  .  Therefore,  the  fall-off  in  intensity  can  be 

compensated  for  by  multiplying  by  ^  ■  Bridou  [22]  has  shovra  that  a  modified  form  of 

can  be  employed  to  correct  for  the  imperfect  periodicity  of  the  oscillations.  This  is 
given  by 

_  _ 


where  m  is  the  order  of  the  interference  and  0c  is  the  critical  angle  as  measured  from  the 
half  maximum  point  of  the  reflectivity  curve. 

The  optical  model  presented  in  this  paper  allows  for  the  reflectivity  of  a  given 
film  to  be  accurately  predicted  from  a  knowledge  of  the  structure.  However,  the  direct 
inversion  of  this  model  to  yield  structure  information  from  the  measured  reflectivity  is 
not  possible  due  to  the  absence  of  phase  information.  Instead,  film  parameters  can  be 
determined  by  fitting  a  multilayer  optical  model  to  the  measured  curve  and  extracting  the 
film  parameters  from  the  model.  This  fit  is  made  by  minimizing  the  mean  square  of  the 
error  between  the  model  and  data,  defined  as 

where  N  is  the  number  of  data  points,  P  is  the  number  of  free  parameters,  and  d  is  the 
experimental  uncertainty  [17].  If  models  with  equal  numbers  of  free  parameters  are  to  be 
used,  the  term  (N-P-1)  can  be  replaced  with  N.  The  minimization  of  this  quantity  is 
effected  by  calculating  the  reflectivity  from  the  model  based  on  an  initial  guess  of  the 
parameters  and  comparing  this  with  the  measured  reflectivity.  A  systematic  search  of 
parameter  space  is  then  performed  according  to  a  particular  algorithm  and  the  reflectivity 
again  calculated  and  compared  with  the  measured  reflectivity.  This  cycle  continues  until 


a  specified  degree  of  convergence  is  reached.  Great  care  must  be  exercised  in  the  design 
of  such  a  routine  as  least-square  fits  to  reflectivity  models  are  known  to  produce  solutions 
which  are  not  unique.  [11]  Moreover,  the  initial  values  chosen  for  the  parameter  search 
may  affect  final  set  of  parameters  which  the  fit  converges  upon.  Because  of  this,  it  has 
been  found  that  the  accuracy  of  a  fit  is  improved  by  providing  initial  values  for  the 
parameters  space  search  which  are  near  to  the  actual  values.  [23] 

Experimental 

To  enable  the  measurement  of  reflectivity  curves  from  growing  films  a 
specialized  system  was  constructed.  A  schematic  representation  of  this  system  is  shown 
in  Figure  13.  The  current  state  of  the  system  is  a  manufactured  chamber  undergoing 
testing,  functioning  acquisition  software,  and  prototyped  analysis  software.  A  second 
x-ray  system,  tailored  to  energy  dispersive  measurement,  is  being  designed  and 
constructed.  It  will  not  initially  be  part  of  the  investigation,  however,  and  will  not  be 
detailed. 

Deposition  Chamber 

The  first  phases  of  this  investigation  will  be  performed  with  a  conventional 
Bragg-Brentano  geometry  diffractometer  having  typical  source-to-sample  distances  on 
the  order  of  tens  of  centimeters.  Because  of  this,  very  particular  restrictions  were  placed 
on  the  design  of  the  deposition  chamber.  Designed  to  use  the  existing  0/20  stage,  the 
chamber  is  small,  rapidly  demountable,  and  reasonably  light-weight.  It  consists  of  two 
independent  substructures  as  shown  in  Figure  14.  The  chamber  itself,  constructed  of  300 
series  stainless  steel,  is  tubular  in  design.  It  consists  of  a  long  flanged  section  in  which 
the  sample  is  housed  and  two  end  caps  which  house  all  of  the  feedthroughs  for  the 
chamber.  Copper-gasketed  conflat  flanges  are  used  at  each  of  these  access  points.  The 
three  sections  of  the  chamber  bolt  together  and  are  sealed  by  Viton  o-rings.  Kapton  film 
is  used  to  seal  the  x-ray  windows  which  run  along  the  sides  of  the  chamber. 


26 


Chamber 


Figure  13.  The  basic  layout  of  chamber,  source  and  detector  showing  Kapton  endows, 
position  of  substrate  and  incident  and  reflected  angles. 
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A  second,  aluminum  construction  attaches  the  chamber  to  the  goniometer  and 
provides  support  and  vibration  isolation  for  the  turbomolecular  pump.  The  chamber  bolts 
to  an  aluminum  plate  which  in  turn  mounts  to  the  goniometer  by  means  of  two  support 
brackets.  The  pump  is  bolted  into  a  mounting  cradle  above  the  chamber.  For  future  use, 
the  deposition  chamber  may  be  adapted  to  another  system  simply  by  redesign  of  this 
mounting  assembly. 

The  sample  stage  rides  on  two  linear  bearing  stages  which  traverse  the  length  of 
the  chamber  and  allow  micrometer  adjustment  of  the  sample  stage  over  a  range  of  two 
inches.  This  travel,  in  addition  to  the  length  of  the  Kapton  windows,  allows  for  a  wide 
range  of  incident  angles  and  substrate  thicknesses.  In  order  to  provide  controlled  growth, 
the  substrate  holder  itself  is  of  copper  block  construction  and  incorporates  temperature 
control  by  means  of  a  Peltier  heating/cooling  element  and  thermocouple  in  addition  to 
water  cooling. 

Deposition  is  presently  performed  by  evaporation,  making  use  of  a  heated 
tungsten  filament  and  thin  foil  sources.  However,  the  chamber  is  designed  to  allow  for 
replacement  of  the  evaporation  source  by  a  small  diameter  sputtering  source  without 
retooling. 

X-ray  System 

For  in  situ  measurement,  an  x-ray  configuration  similar  to  that  used  in  the  current 
system  for  energy  dispersive  measurements  is  anticipated.  This  configuration  is 
described  by  Windover,  et.  al.  [1] 

[We  used  a]  20  kV,  20  mA  (400  Watt)  0.4  mm  x  12  mm  copper  anode  sealed  source  X-ray 
tube  on  a  commercial  Scintag  X-ray  dif&actometer.  The  source-to-sample  and  detector-to 
sample  distance  was  286  mm.  Source  and  scatter  0.05  mm  collimator  slits  were  used  on  the 
source  side,  providing  0.043°  beam  divergence  in  the  reflection  plane.  One  detector  0.05 
mm  collimator  slit  was  used  providing  a  fixed  angle  with  0.01°  width.  A  nominal  600-pm 
thick,  Si  (100)  single  side  polished  wafer  was  used  in  the  primary  beam  path  to  filter  the  Cu 
ka  and  kp  lines  firom  the  incident  beam  spectruia 

Through  the  use  of  the  silicon  wafer  as  a  filter,  the  intensity  of  the  ka  and  kP  peaks 
relative  to  the  broad  spectrum  brehmstrahlung  is  reduced.  This  enables  higher  beam 
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intensity  and  increased  count  rates  across  the  spectrum  without  detector  saturation.  The 
more  even  spectrum  also  improves  the  accuracy  of  the  intensity  normalization. 

Figure  1 1  shows  the  spectrum  of  the  incident  beam  and  the  impact  of  normailzation  on 
the  reflectivity  data. 

Previous  work  with  flxed-angle,  energy  dispersive  X-ray  reflectivity  has  involved 
using  a  high  voltage  source.  The  approach  adopted  for  this  work  will  be  to  use  lower 
voltage  sources,  copper  or  chromium  operating  at  20  kV  to  produce  only  the  region  of  the 
reflectivity  spectrum  needed  to  recover  the  thickness  or  density  information. 
Semiconductors  are  prone  to  device  degradation  by  X-ray  exposure  due  to  the  creation  of 
charge  traps  caused  by  X-ray  generated  defects  [1].  By  using  a  lower  power  X-ray 
source  and  minimizing  the  exposure  time,  this  measurement  method  can  be  made  suitable 
for  routine  semiconductor  device  characterization. 


Data  Acquisition 


Reflected  intensity  is  measured  by  a  Peltier  cooled,  Silicon  Kevex  detector.  This 
device  has  an  energy  resolution  of  180  eV  FWHM  at  5.6  keV,  which  is  the  equivalent  of 
a  scan  resolution  of  a  few  hundredths  of  a  degree.  The  relation  between  the  two 
resolutions  is  determined  fi-om  the  resulting  resolutions  in  . 
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A  Nucleus  PC  A II  multichannel  analyzer  set  to  256-channel  mode  is  currently 
used  for  captunng  spectral  data.  For  the  range  of  energies  typically  considered,  this 
allows  for  a  resolution  in  greater  than  that  imposed  by  the  detector.  In  previous 
energy  dispersive  measurements  the  linearity  and  calibration  was  set  using  7  data  points; 
Cr  ka  and  k(3,  Cu  ka  and  kp,  and  Mo  ka  and  kp,  and  the  20  kV  edge  in  the  data. 

Data  from  the  card  is  stored  to  a  file  using  a  routine  written  in  C.  (See  q)pendix 

1) 
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Computational  Analysis 


After  a  sufficient  number  of  counts  have  accumulated  on  the  MCA,  the  value 
from  the  card  is  read  to  a  file.  This  data  is  then  transferred  into  an  analysis  routine.  For 
rapid  monitoring  of  thickness,  a  Fourier  transform  may  be  used.  This  is  implemented 
within  Lab  View  using  a  windowed  FFT  and  a  peak-fitting  function.  The  pretransform 
and  normalization  can  both  be  performed  directly  within  Lab  View. 

For  more  in-depth  study,  a  fitting  routine  based  on  the  optical  transfer  matrix  is 
utilized.  (The  most  current  version  of  this  software  and  is  detailed  in  appendix  2.)  In  this 
routine,  each  successive  data  set  is  represented  as  another  layer  in  the  model,  having 
unique  parameters.  A  new  multilayer  model  having  a  number  of  layers  equal  to  the 
number  of  data  sets  which  have  been  recorded  is  calculated  for  each  data  set.  The 
parameters  of  the  new  layer  are  determined  from  a  fit  of  this  model  to  the  experimental 
data.  Multilayer  fits  and  the  challenges  associated  with  them  have  been  discussed.  Of 
particular  concern  is  the  convergence  of  the  fit  to  inaccurate  results.  However, 
developing  the  model  by  sequential  fits  made  during  the  film  growth  minimizes  the  non¬ 
uniqueness  problem  by  reducing  the  number  of  fi’ee  parameters.  Only  six  parameters  are 
required  for  each  fit:  thickness,  density,  sur&ce  roughness,  interface  roughness,  5  and  p. 
Further,  the  chance  of  ftilse  convergence  is  minimized  by  taking  initial  values  for  the 
parameter  search  fi-om  the  fitted  values  of  the  previous  layer.  A  flow  chart  of  this  process 
is  shown  in  Figure  15. 

The  Levenberg-Marquardt  algorithm  is  used  to  fit  the  data.  Reflectivity  data  does 
not  have  a  linear  dependence  on  its  parameters  as  given  by 

Therefore,  a  general  nonlinear  fitting  algorithm  is  used.  The  Levenberg-Marquardt  has 
the  advantage  of  being  able  to  converge  on  the  correct  minimum  from  far  away  while 
also  being  able  to  rapidly  converge  to  the  answer  when  the  initial  guess  is  close  to  the 
actual  value.  [24]  The  code  for  this  procedure  is  borrowed  largely  from  Press,  et.  al. 

[25]. 
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Figure  15.  A  flow  chart  illustrating  the  operation  of  the  measurement  software. 
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Summary 


Through  theoretical  exposition,  the  feasibility  of  energy  dispersive  x-ray 
reflectivity  for  the  measurement  of  growing  thin  films  has  been  shown.  A  model  for  the 
reflectivity  based  on  Fresnel  theory  has  been  detailed.  This  model  accounts  for  the 
effects  of  the  thickness,  roughness,  and  density  of  layers  and  multilayers  on  the 
reflectivity  of  a  sample.  Methods  for  the  extraction  of  film  parameters  from  the 
reflectivity  based  on  optical  and  scattering  formalisms  have  been  detailed.  In  order  to 
study  film  growth  by  energy  dispersive  reflectivity  method  a  system  has  been  designed 
and  constructed.  Software  has  been  written  in  Lab  View  and  C-h-  to  acquire  and  analyze 
data  using  both  the  Fourier  transform  method  and  the  least-squares  fit  to  the  optical 
model.  Preliminary  testing  of  this  system  is  now  underway. 
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Appendix  1:  C++  code  to  store  the  data  on  the  PCA II  card  to  a  file 


/*  PROGRAM:  Cdemo.c  >  Cdemo.exe  */ 

/*  COMPILER:  Borland  Turbo  C  ver  2.0  */ 

/*  DATE:  Fall  1988  */ 

/*  COMPANY:  Nucleus  Inc.  */ 

/♦PURPOSE:  read  PCA  Board  RAM,  */ 

/*  then  and  write  an  ascii  */ 

/♦  disk  file  called  ram.dat.  */ 

/*  HARDWARE:  PCA-H  */ 

/*  */ 

/*  modified  for  MS  visual  C++  */ 

/*  June  30,  1999  */ 

/*  Jason  Summers  */ 

/* - Supporting  Include  Files - */ 

#include  <process.h> 

#include  <memory.h> 

#include  <stdio.h> 

#include  <windows.h> 

^include  <conio.h> 

/* - Supporting  Include  Files - */ 

/* - Macros - - */ 

#define  TRUE  1 
#define  FALSE  0 
#define  NIL  0 

#define  NUMB_CHAN  1024  /*  number  of  channels  to  process  */ 

#define  PORTNUMB  0x0  leO  /*  for  default  dip  switch  setting  */ 

^define  DATAFILENAME  "ram.dat"  /*  name  of  data  (disk)  file  to  write  to  */ 
#define  SOURCEDESC  2  /*  Global  Descriptor  Table  for:  source  24  bit  address 

*/ 

#define  TARGETDESC  3  /*  Global  Descriptor  Table  for:  target  24  bit  address 

*/ 

#define  begfimct  { 

#define  endfunct  } 

#defme  begfor  { 

#define  endfor  } 

#define  begwhile  { 

#define  endwhile  } 

#defme  begif  { 

#define  endif  } 
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#define  begelse  { 

#defme  endelse  } 

#define  begelseif  { 

#define  endelseif } 

#defme  begdo  { 

#define  enddo  } 

#defme  begswitch  { 

#define  endswitch } 

#define  BOOL  unsigned  char 

#define  SHORT  int 

#define  USHORT  unsigned  short 

#define  LONG  long 

#define  ULONG  unsigned  long 

#define  CHAR  char 

#defineUCHAR  unsigned  char 

/* - Macros - */ 

/* .  - . - - - DATA  TYPES 

— — . .  */ 

/* - CONTROL  STRUCT - */ 

typedef  struct 

{ 

USHORT  cardport;  /*  number  of  lowest  used  PCA-II  I/O  port  */ 

USHORT  cardtopport;  /*  highest  port  address  used  by  PCA-II  (always: 


cardport  +  2)  */ 

}  CONTROL; 

/* - CONTROL  STRUCT - */ 

/* - GDT  STRUCT - */ 

typedef  struct  { 

unsigned  int  seg_limit;  /*  segment  size  limit  of  move  */ 


unsigned  int  base_lo_word;  /*  lower  16  bits  of  24  bit  address  */ 

unsigned  char  base_hi_byte;  /*  upper  8  bits  of  24  bit  address  */ 

unsigned  char  data_access_rights;  /*  access  rights  =  0x93  under  DOS  3.0  */ 
unsigned  int  data_reserved;  /*  reserved  for  use  by  interrupt  */ 

}  GDT; 

/* - GDT  STRUCT - */ 

/*  - - =  DATA  TYPES 

/* - Global  Variables - */ 

ULONG  laiTay[NUMB_CHAN];  /*  count  array  */ 

UCHAR  roi[NUMB_CHAN];  /*roi  array  */ 

CONTROL  control;  /*  data  structure  */ 

/* - *! 
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j:if.  GLOBAL  DECRIP  TOR  TABLE 

^f:*****^  *3|c:|c:tct)fc)|(9(c%%)(c3fc3|c;4c3|c  if  j 

GDTgclt[]={ 

{  0,  0,  0,  0,  0  },  /*  dummy  table  for  interrupt  */ 

{  0,  0,  0,  0,  0  ),  /*  dummy  table  for  interrupt  */ 

{  0,  0,  0,0x93,  0  },  /*  table  for  the  source  of  move  */ 
{  0,  0,  0,0x93,  0  },  /*  table  for  the  target  of  move  */ 

{  0,  0,  0,  0,  0  ),  /*  dummy  table  for  interrupt  */ 

{  0,  0,  0,  0,  0  )  /*  dummy  table  for  interrupt  */ 

}; 

/*  ^**t******************  global  DECRIPTOR  TABLE 


ic  }|c  ^  t  ***  t  sK  t  9ft  It  jft  ^ 

/* - Global  Variables - */ 

/* - Function  Prototypes - */ 


void  WINAPI  putportbyte(register  SHORT  address,  UCHAR  value); 
UCHAR  WINAPI  getportbyte(register  SHORT  address); 

void  get_data_PCAII(void); 


/* - Function  Prototypes - */ 

/* - Main  Program  Start: - */ 

void  main(void) 

{ 

register  SHORT  i;  /*  channel  loop  index  */ 

FILE  *fp;  /*  file  pointer  */ 

CHAR  ch;  /*  temp  char  */ 

SHORT  ret  =  0;  /*  var  for  function  return  code  */ 


memset(  &larray[0],  NIL,  (NUMB_CHAN  *  sizeof(ULONG)) );  /*  init  the  long  int 
array  to  zero  */ 


control.cardport  =  PORTNUMB;  /*  save  for  CONTROL  */ 
control.cardtopport  =  PORTNUMB  +  2;  /*  save  for  CONTROL  */ 

ch  =  getportbyte(O);  /*  get  data  on  card,  (first  byte)  */ 
ret  =  ch;  /*  save  so  we  can  restore  later  on  */ 

ch-H-;  /*  add  one  to  the  read  value  */ 

putportbyte(0,ch);  /*  write  out  the  incremented  data  */ 

if  (getportb5^e(0)  !=  ch)  /*  PCA-II  card  not  present  */ 
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begif 

puts("No  PCA-II  Card  was  Found  at  port  OleOh."); 
endif 

else  /*  a  PCA-II  Card  was  present,  so  use  PCA-II  card  */ 
begelse 

putportbyte(0,  (UCHAR)  ret);  /*  restore  the  changed  data  */ 
get_data_PCAII();  /*  read  data  from  the  PCA-II  card  */ 
endelse 


/*  NOTE:  Now  we  have  both  channel  counts  and  ROI  data  copied  */ 

/*  into  easy  to  use  arrays.  So  now  you  can  do  your  */ 

/*  calculations  on  the  data  in  the  card.  */ 

/*  channel  counts  in:  laiTay[]  (a  long  int  array)  */ 

/*  ROI  data  in:  roi[]  (a  char  array)  */ 

/*  NOTE:  (if  an  ROI  byte  is  not  zero  then  that  channel  is  part  of  an  ROI  */ 

fp  =  fopen(DATA_FILE_NAME,"w");  /*  note:  data  file  name,  write  mode  */ 

if  (fp  !=  NULL)  /*  there  was  no  problem  opening  the  disk  file  */ 
begif 

fprintf(fp,"Channel  Counts  ROI  \n");  /*  some  header  text  */ 
fprintf(fp," - \n"); 

/* - Output  Loop  to  Disk  File - */ 

for  (i  =  0;  i  <  NUMB_CHAN;  i++) 
begfor 

if  (  (i  %  512)  =  0)  printf("writing  record:  %d\n",i); 

ret  =  fprintf(fp,"%4d  %81d  %4d\n",i,larray[i],roi[i]); 

if  (ret  =  -1)  /*  test  for  i/o  (disk  full)  error  */ 
begif 

puts("ERROR:  Disk  error.  Press  any  key  to  continue."); 


getcharO;  /*  pause  for  a  key  */ 
exit(-9);  /*  exit  to  DOS  */ 

endif 

endfor 

/* - Output  Loop  to  Disk  File - */ 


i--;  /*  adjust  loop  index  */ 
printf("writing  record:  %d\n",i); 
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fclose(^);  /*  close  the  disk  file  */ 

printf("Write  of  ascii  file:  %s  complete.\n",DATA_FILE_NAME); 
endif 
else 
begelse 

putsC'ERROR;  Opening  Disk  File.'');  /*  print  error  message  */ 
endelse 

endfiinct  /*  end  of  main  program  */ 

/* - Main  Program  End: - */ 


/*  needed  subroutines  ('functions'  in  C  Language)  */ 

/*  This  function:  writes  data  to  PCA-II  I/O  mapped  MCA  boards  */ 

/*  =========_========_==============_^ 

*/ 

void  WINAPI  putportbyte(SHORT  address,  UCHAR  value) 
begfimct 

USHORT  port  =  control,  cardtopport;  /*  port  number  plus  two  */ 

_outp(port--,  address  );  /*  low  address  byte  */ 

address  =  (address  »  8);  /*  move  the  hi  b3^e  into  the  lo  position  */ 

_outp(port— ,  address  );  /*  high  address  byte  */ 

_outp(port,  value );  /*  put  the  data  byte  */ 

endfunct 
/* 


*/ 

/*  This  function:  reads  data  from  PCA-II  I/O  mapped  MCA  boards  */ 
/*  ==— =====^=====_== 
*! 

UCHAR  static  WINAPI  getportbyte(SHORT  address) 
begfunct 

USHORT  port  =  control.cardtopport;  /*  port  number  plus  two  */ 

_outp(port~,  (UCHAR)  address);  /*  low  address  byte  */ 
_outp(port-,  address  »  8);  /*  high  address  byte  */ 

retum(  _inp(port) );  /*  get  the  data  byte  */ 

endfunct 
/* 


*/ 
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/*  This  function;  gets  count  and  ROI  data  from  a  PCA-II  card  */ 

/*  - - - 

*/ 

void  get_data_PCAII(void) 
begfunct 

register  SHORT  i;  /*  loop  index  */ 

register  SHORT  cardbyte;  /*  byte  loop  index  */ 

CHAR  *byteptr;  /*  pointer  to  count  buffer  */ 

cardbyte  =  0;  /*  init  RAM  byte  index  */ 

/* - Loop  to  fill  the  long  integer  array - */ 

for  (i  =  0;  i  <  NUMB_CHAN;  i++) 
begfor 

byteptr  =  (UCHAR  /*far*/  *)  (&larray[i]);  /*  set  byte  ptr  to  the  long  int  array  element 

*/ 

*byteptr++  =  getportbyte(cardbyte-H-);  /*  get  low  byte  of  24  bit  channel  value  */ 
*byteptr-H-  =  getportbyte(cardbyte++);  /*  get  mid  byte  of  24  bit  channel  value  */ 
*byteptr  =  getportbyte(cardbyte-H-);  /*  get  hi  byte  of  24  bit  channel  value  */ 
roi[i]  =  getportbyte(cardbyte-H-);  /*  get  ROI  byte  */ 
endfor 

/* - Loop  to  fill  the  long  integer  array - */ 

endfunct 

/* 


Appendix  2:  C++  code  which  makes  a  sequential  fit  of  a  growing  film  to  a  multilayer 
model 


/*  Program  to  to  test  the  curve  fitting  routine;  testfit.cpp  */ 

#include  <stdio.h> 

^include  <math.h> 

#include  <stdlib.h> 

typedef  struct  FCOMPLEX  (float  r;  float  i;}  fcomplex; 

#define  NPT  256  /*  number  of  channels  on  MCA*/ 
#define  MA  6  /*  number  of  parameters  per  layer  */ 

#define  SPREAD  0.001 
#define  SQR(a)  ((a)*(a)) 

#define  SW^(a,b)  (float  temp=(a);(a)=(b);  (b)=temp;} 


/* - Function  Prototypes - */ 

void  covsrt(float  **covar,  int  ma,  int  *Iista,  int  mfit); 

void  gaussj  (float  **a,  int  n,  float  **b,  int  m); 

void  mrqcof(float  *x,  float  *y,  float  *sig,  int  ndata,  float  *a,  int  ma, 
int  *lista,  int  mfit,  float  **alpha,  float  *beta, 
float  *chisq,  void(*refcurve)  (float, float*, float,float*, 
float,  int)); 

void  mrqmin(float  *x,  float  *y,  float  *sig,  int  ndata,  float  *a,  int  ma, 
int  *lista,  int  mfit,  float  **covar,  float  **alpha, 
float  *chisq,  void(*refcurve)  (float,float*,float,float*, 
float,  int),float  *alamda); 

void  refcurve(float  x,  float  *a,  float  y,  float  *dyda,  float  theta,  int  N); 

void  refmatrix(int  N,  int  layer,  float  d,  float  lambda,  float  delta  1,  float  delta2, 

float  thetal,  float  theta2,  fcomplex  oldmatrix[l][l], 
fcomplex  newmatrix[l][l]); 

fcomplex  Complex(float  re,  float  im); 

float  Cabs(fcomplex  z); 

fcomplex  Conjug(fcomplex  z); 

fcomplex  Csub(fcomplex  a,  fcomplex  b); 

fcomplex  Cmul(fcomplex  a,  fcomplex  b); 
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fcomplex  Cdiv(fcomplex  a,  fcomplex  b); 
fcomplex  Cadd(fcomplex  a,  fcomplex  b); 
fcomplex  RCmul(float  x,  fcomplex  a); 

float  *vector(int  nl,  int  nh); 

int  *ivector(int  nl,  int  nh); 

void  free_vector(float  *v,  int  nl,  int  nh); 

void  free_ivector(int  *v,  int  nl,  int  nh); 

void  nrerror(char  *eiTor_text); 

float  **matrix(int  nrl,  int  nrh,  int  ncl,  int  nch); 

void  free_matrix(float  **m,  int  nrl,  int  nrh,  int  ncl,  int  nch); 


/♦ 


Function  Prototypes - */ 


/* - Global  Variables - */ 

float  theta=l; 
intN=3; 
fcomplex  one; 
one.r  =  1;  one.i=0; 

/* - Global  Variables - */ 


int  main(void) 

{ 

FILE  ^normalization,  *data; 

int  i,  idum=(-911),itst,k,mfit,*lista,  bin[NPT+l],  bin2[NPT+l]; 
unsigned  int  larray[NPT+l],  normarray[NPT+l]; 
float  alamda,chisq,ochisq,  *x,  *y,  *  sig,  *  *covar, 

**  alpha; 

static  float  a[MA+l]= 

{0.0,5.0,2.0,3.0,2.0,5.0,3.0};  /*  a[l]=theta  in;  a[2]=lambda; 

a[3]=delta; 

a[4]=beta;  a[5]=thickness; 

a[6]=top  roughness  */ 

static  float  gues[MA+l]= 

{0.0,4.5,2.2,2.8,2.5,4.9,2.8}; 


/*  read  the  normalization  file  in  */ 
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fopen("c:\\norm.dat","r"); 

i-0; 

while  (i<NPT) 

{ 

fscanf(normali2ation,  ''%i  %u\n'',&bin[i],&normarray[i]); 

} 

fclose  (normalization); 

/*  read  the  data  file  in  */ 
fopen("c:\\ram.dat","r"); 
i=0; 

while  (i<NPT) 

{ 

fscanf(data,"%i%u\n",&bin2[i],&larray[i]); 

} 

fclose  (normalization); 


lista=ivector(  1  ,MA); 
x=wector(  1  ,I^T); 
y=vector(  1  ,NPT); 
sig=vector(  1  ,NPT); 
covar=matrix(  1  ,MA,  1  ,MA); 
alpha=matrix(  1  ,MA,  1  ,MA); 
for  (i=l;i<=NPT;i++)  { 

x[i]=24*i+34;  //set  energy  scale 

y[i]=larray[i]/normarray[i];  //perform  normalization 
sig[i]=SPREAD*y[i]; 

} 

mfit=6; 

for  (i=l;i<=mfit;i-H-)  lista[i]=l; 
alamda  =  -l; 

for  (i=l;i<=MA;i-H-)  a[i]=gues[i]; 

mrqmin(x,y,sig,NPT,a,MA,lista,mfit,covar,  alpha,  &chisq,  refcurve,  &alamda); 
k=l; 
itst=0; 
while  (itst  <  2) 

{ 

printf(7n"); 

k-H-; 

ochisq=chisq; 

mrqmin(x,y,sig,NPT,a,MA,lista,mfit,covar,alpha,&chisq,refcurve,&alamda); 

if  (chisq  >  ochisq)  itst=0; 

else  if  (fabs(ochisq-chisq)  <  0. 1)  itst-H-; 
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alamda=0.0; 

mrqmin(x,y,sig,NPT,a,MA,lista,mfit,covar,alpha,&chisq,refcurve,&alamda); 

printf("\nUncertainties:\n"); 

for  (i=l;i<=6;i-H-)  printf)["%9.4f',sqrt(covar[i]  [i])); 

printf(''\n"); 

ffee_matrix(alpha,  1,MA,  1,MA); 
free_matrix(covar,  1  ,MA,  1  ,M  A); 
free_vector(sig,  1  ,NPT); 
free_vector(y,  1  ,NPT); 
ff  ee_vector(x,  1  ,NPT); 
free_ivector(Iista,  1  ,MA); 
return  (0); 

} 


/*  Levenberg-Marquardt  minimization  of  chi-squared  for  a  fit  to  a 
non-linear  function  with  parameters  a[l..ma]  */ 

*! 

void  mrqmin(float  *x,  float  *y,  float  *sig,  int  ndata,  float  *a,  int  ma, 
int  *lista,  int  mfit,  float  **covar,  float  **alpha, 
float  *chisq,  void(*refcurve)  (float, float*, float*, float*, 
float,  int),float  *alamda) 


int  k,kk,j,ihit; 

static  float  *da,*atry,**oneda,*beta,ochisq; 


if  (*alamda  <  0.0)  { 

oneda=matrix(  1 ,  mfit,  1,1); 
atry=vector(  1  ,ma); 
da=vector(l,ma); 
beta=vector(l,ma); 
kk=mfit+l; 

/*  check  to  see  if  the  lista  values  are  reasonable  */ 


forO=l;j<=ma;j-H-)  { 
ihit=0; 

for  (k=  1  ;k<=mfit  ;k-i-+) 

if  (lista[k]  =  j)  ihit-H-; 
if  (ihit^^^O) 

lista[kk-H-]=j; 
else  if  (ihit  >1) 
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nrerror  ("Bad  LISTA  permutation  in  MRQMIN-l"); 


} 

if  (kk  !=  (ma+1))  nrerror  ("Bad  LISTA  permutation  in  MRQMIN-2"); 
/*  okay  so  proceed  */ 

*alamda=0.001; 

mrqcof(x,y,sig,ndata,a,ma,lista,mfit,alpha,beta,chisq,refcurve); 

ochisq=(*chisq); 

} 

for(j=l;j<=mfit;j++)  { 

for  (k=l;k<=mfit;k-H-) 

covar  0][k]=alpha0][k]; 
covar[j][j]=alpha[j]  D]*(  1  •0+(*alamda)); 
oneda[j][l]=beta[j]; 

} 

gaussj(covar,mfit,oneda,  1); 

for  (j=l;j<=mfit;j-H-)  da[j]=oneda[j][l]; 

if  (*alamda  —  0.0) 

{ 

covsrt(covar,ma,lista,mfit); 
free_vector(beta,  l,ma); 
free_vector(da,  l,ma); 
ff ee_vector(atry,  1 ,  ma); 
free_matrix(oneda,  1  ,mfit,  1,1); 
return; 

} 

for  (j=l;j<=ma;j-H-)  atry(j]=a|j]; 
for  (j=l;j<=mfit;j-H-)  atry[lista[i]]=a[listaO]]+da[j]; 
mrqcof(x,y,sig,ndata,atry,ma,lista,mfit,covar,da,chisq,refcurve); 
if(*chisq  <  ochisq)  { 

*alamda  *=0.1; 
ochisq=(*chisq); 
for  (j=l;j<=nifit;j++)  { 

for  (k=l;k<=mfit;k-i-+) 

alpha[j]  [k]=covar[j]  [k]; 
beta[j]=da|j]; 
a[lista[j]]=atry[lista[j]]; 

} 

}  else  { 

*alamda  *=  10.0; 

*chisq=ochisq; 

} 

return; 

} 
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/*  This  performs  math  duties  for  the  routine  mrqmin  */ 

void  mrqcof  (float  *x,  float  *y,  float  *sig,  int  ndata,  float  *a,  int  ma, 
int  *lista,  int  mfit,  float  **alpha,  float  *beta, 
float  *chisq,  void(*refcurve)  (float*, float*,float*, float*, 
float,  int)) 


int  k,j,i; 

float  ymod,wt,sig2i,dy,*dyda; 

dyda=vector(  l,ma); 
for(j=l;j<=mfit;j++){ 

for  (k=l;k<=j;k-H-)  alphalj][k]=0.0; 
beta[j]=0.0; 

} 

*chisq=0.0; 

for  (i=l;i<=ndata;i++)  { 

//(*refcurve)(x[i],a,&ymod,dyda,theta,N); 

sig2i=1.0/(sig[i]*sig[i]); 

dy=y[i]-ymod; 

for  (j=l;j<=mfit;j-H-)  { 

\vt=dyda[lista(j]]*sig2i; 
for  (k=l;k<=j;k++) 

alpha[j][k]  +=  'wt*dyda[lista[k]]; 
beta|j]  +=  dy*wt; 

} 

(*chisq)  +=  dy*dy*sig2i; 

} 

for0=2;j<=mfit;j++) 

for  (k=  1  ;k<=j- 1  ;k-H-) 

alpha[k][j]  =  alphalj][k]; 
free_vector(dyda,  1 ,  ma) ; 


/*  This  sorts  covariance  matrix  */ 

j* - - - - - - 

void  covsrt(float  **covar,  int  ma,  int  *lista,  int  mfit) 

{ 


.*/ 


int  i,j; 
float  swap; 
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for  (j=l;j<nia;j++) 

for  (i=j+l;i<=ma;i++)  covar[i][j]=0.0; 
for  (i=l;i<mfit;i++) 

for(j=i+l;j<=mfit;j++)  { 

if  (listafj]  >  lista[i]) 

covar[lista|J]]  [lista[i]]= 
covar[i][j]; 

else 

covar[lista[i]][Iista|j]]= 

covar[i]Ij]; 

} 

swap=covar[l][l]; 

for(j=l;j<=ma;j++){ 

covar  [1]0]  =  covar  [j][j]; 
covar[j][j]=0.0; 

} 

covar[lista[  1  ]][lista[  1  ]]=swap; 
for  (i=2;j<=mfit;j-H-)  covar[lista|j]] 
[lista[j]]=covar[l]lj]; 
for  (j=2;j<=ma;j-H-) 

for  covar[i][j]= 

covar[j][i]; 


} 


/*  Gauss-Jordon  sol'n  to  lin.  eq.  */ 

I* _ _ _ */ 

void  gaussj  (float  **a,  int  n,  float  **b,  int  m) 


int  *indxc,*indxr,*ipiv; 
int  i,icol,irow,j,k,l,ll; 
float  big,dum,pivinv; 

indxc=ivector(l,n); 
indxr=ivector(  1  ,n); 
ipiv=ivector(  1 ,  n); 
for  0=l;j<=n;j+-t-)  ipiv[j]=0; 

for  (i=l;i<=n;i-H-)  ( 

big=0.0; 

for 

if  (ipivU]  !=  1) 
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for  (k=l;k;<=n;k-H-)  { 

if(ipiv[k]  =  0){ 
if  (fabs(a[j][k])  >=  big)  { 

big=fabs(a[j][k]); 

irow=j; 

icol=k; 

} 

}else  if  (ipiv[k]  >  1)  nrerror  ("GAUSSJ;  singular 

matrix-1 "); 

} 

-H-(ipiv[icol]); 

if  (irow  !=  icol)  { 

for  (l=l;l<=n;l++) 

SWAP  (a[irow][l],a[icol][l]) 
for  (l=l;I<=m;l++) 

SWAP  ^[irow][l],b[icol][l]) 

} 

indxr[i]=irow; 
indxc[i]=icol; 
if  (a[icol][icol]  =  0.0) 

nrerror  ("GAUSSJ:  singular  matrix-2"); 
pivinv=l  .0/a[icol][icol]; 
a[icol][icol]=1.0; 

for  (l=l;l<=n;l-H-)  a[icol][l]  *=pivinv; 
for  (l=l;l<=m;l++)  b[icol][l]  *=pivinv; 
for(ll=l  ;ll<=n;ll++) 
if  (11  !=  icol)  { 

dum=a[ll][icol]; 

a[ll][icol]=0.0; 
for(l=l;l<=n;l++) 
a[ll][l]  -=  a[icol][l]*dum; 
for  (l=l;l<=m;l-H-) 
b[ll][l]  -=  b[icol][l]*dum; 

} 

} 

for(l=n;l>=l;l-)  { 

if(indxr[l]  !=indxc[l]) 

SWAP(a[k][indxr[l]],a[k][indxc[l]]); 

} 

ff ee_ivector(ipiv,  1  ,n); 
ff  ee_ivector(indxr,  1  ,n); 
ff  ee_ivector(indxc,  1  ,n) ; 
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/*  This  calculates  reflectivity  for  N  layers  */ 
/* 


void  refcurve  (float  x,float  *a,float  y,float  *dyda,float  theta, int  N) 

{ 

intj; 

float  *deltaA,  yplus; 
dyda=vector(  1  ,NPT); 
deltaA=vector(  1  ,MA); 


fcomplex  oldmatrix[l][l],  newmatrix[l][l],  r,  R; 
float  d, lambda,  deltal,  delta2,  thetal,  theta2; 


oldmatrix  [0][0]  =  Complex(l,0); 
oldmatrix  [1][1]  =  Complex(l,0); 
oldmatrix  [1][0]  =  Complex(0,0); 
oldmatrix  [0][1]  =  Complex(0,0); 


j=i; 

for  (j=l;j<=N;j4-f) 

{ 

layer=j; 

refmatrix  (N,layer,d,  lambda,deltal,delta2,thetal,theta2,oldmatrix,newmatrix); 

} 

r=Cdiv(newmatrix[0][  1  ],newmatrix[  1  ]  [  1  ]); 
y=Gmul(r,r);  /*  the  reflectivity  */ 


/*  the  following  takes  a  numerical  derivative  of  theoretical 
y(lambda)  w.r.t  a  given  parameter*/ 

for  (j=l;j<=MA;j-H-)  deltaA[j]=  00001  *a[j]; 


for  (j=l;j<=MA;j++) 

{ 

if  0=1) 

{ 

forO=l;j<=Ny++) 

{ 
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layer=j; 

refmatrix  (N, layer, d,  lambda, deltal ,delta2,thetal ,theta2,oldmatrix,newmatnx); 

} 

r=Cdiv(newmatrix[0][l],newmatrix[l][l]); 

y=Cmul(r,r); 

} 

aO]=a[j]+deltaAO]; 
for  0=i;j<=N;j++) 

{ 

layer=j; 

refmatrix  (N,layer,d,  lambda,deltal,delta2,thetal,theta2,oldmatrix,newmatnx); 

} 

r=Cdiv(newmatrix[0]  [  1  ]  ,newmatrix[  1  ]  [  1  ]); 
yplus=Cmul(r,r); 

aIj]=a[j]-deltaA[j]; 

for 

{ 

layer=j; 

refmatrix  (N,layer,d,  lambda,deItal,delta2,thetal,theta2,oldmatnx,newmatnx); 

} 

r=Cdiv(newmatrix[0]  [  1  ]  ,newmatrix[  1  ]  [  1  ]); 
y=Cmul(r,r); 

dyda|j]=((yplus-y)/deItaAD]); 

} 


) 


/*  Calculates  the  optical  matrix  for  refcurve  */ 

/* 

=*l 

void  refmatrix  (int  N,  int  layer,  float  d,  float  lambda,  float  deltal, 

float  delta2,float  thetal,  float  theta2, 

fcomplex  oldmatrix[l][l],fcomplex  newmatrix[l][l]) 

{ 

fcomplex  phi,  pi,  p2,  t,  ttilda,  r,  rtilda,  nl,  n2,  A[l][l]; 
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float  betal,  beta2; 


nl.r=  1-deltal;  n2.r=  l-delta2;  nl.i=-betal;  n2.i=-beta2; 
theta2  =  acos  ((nl.r/n2.r)*cos(thetal)); 

p  1  =RCmul(sin(theta  1 ),  n  1 ); 

p2=RCmul(sin(theta2),n2); 

r=Cdiv(Csub(p  1  ,p2),Cadd(p  1  ,p2)); 

t=Cdiv(RCmul(2,p  1  ),Cadd(p  1  ,p2)); 

rtilda=RCmul(- 1  ,r); 

ttilda=Cdiv(Cmul(t,p2),p  1 ); 

phi=RCmul((2*3.1415*d*sin(thetal))/lambda,nl); 


if  (layer!=N) 

{ 

A[0][0]=Cmul(Csub(Cmul(t,ttilda),Cmul(r,rtilda)),exp(2*i*phi)); 

A[0][l]=r; 

A[l][0]=RCmul(-l,rtilda); 

A[l][l]=l; 

} 

else 

{ 

A[0]  [0]=Cdiv(one,t); 

A[0][l]=Cdiv(r,t); 

A[l][0]=Cdiv(r,t); 

A[l][l]=Cdiv(one,t); 

} 

newmatrix[0][0]=Cadd(Cmul(oldmatrix[0][0],A[0][0]),Cmul(oldmatrix[0][l],A[l 

][0])); 

newmatrix[0][  1  ]=Cadd(Cmul(oldmatrix[0]  [0],  A[0][  1  ]),Cmul(oldmatrix[0][  1  ],  A[  1 

][!])); 

newmatrix[  1  ]  [0]=Cadd(CmuI(oldmatrix[  1  ]  [0],  A[0][0]),Cmul(oldmatrix[  1  ][  1  ],  A[  1 

][0])); 

newmatrix[  1  ]  [  1  ]=Cadd(Cmul(oldmatrix[  1  ]  [0] ,  A[0]  [  1  ]),Cmul(oldmatrix[  1  ]  [  1  ],  A[  1 

][!])); 


} 


/*= 

*/ 
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/*  The  following  allows  C-H-  to  handle  complex  numbers  */ 
/*  ■■ 

_*/ 


fcomplex  Complex(float  re,float  im) 

{ 

fcomplex  c; 
c.r=re; 
c.i=im; 
return  c; 

} 


float  Cabs(fcomplex  z) 

{  float  x,y,ans,temp; 
x=fabs(z.r); 
y=fabs(z.i); 
if(x=0.0) 

ans=y; 

else  if  (y=0.0) 
ans=x; 
else  if  (x>y)  { 

temp=y/x; 

ans=x*sqrt(l  .0+temp*temp); 

}else{ 

temp=x/y; 

ans=y*sqrt(l  .0+temp*temp); 

} 

return  ans;} 

fcomplex  Conjug  (fcomplex  z) 

{ 

fcomplex  c; 
c.r=z.r; 
c.i=-z.i; 
return  c; 

} 

fcomplex  Csub(fcomplex  a,fcomplex  b) 

{ 
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fcomplex  c; 
c.r=a.r-b.r; 
c.i=a.i-b.i; 
return  c; 

) 

fcomplex  Cmul(fcomplex  a,  fcomplex  b) 

{ 

fcomplex  c; 
c.r=a.r*b.r-a.i*b.i; 
c.i=a.i*b.r+a.i*b.i; 
return  c; 

} 

fcomplex  Cdiv(fcomplex  a,fcomplex  b) 

{ 

fcomplex  c; 
float  r,den; 

if  (fabs(b.r)  >=  fabs(b.i))  { 
p=b.i/b.r; 
den=b.r+r*b.i; 
c.  r=(a.  r+r*a.  i)/den; 
c.i=(a.i-r*a.r)/den; 

}else( 

r=b.r/b.i; 

den=b.i+r*b.r; 

c.r=(a.r*r+a.i)/den; 

c.i=(a.i*r-a.r)/den; 

} 

return  c; 


fcomplex  Cadd(fcomplex  a,fcomplex  b) 


{ 

fcomplex  c; 
c.r=a.r+b.r; 
c.i=a.i+c.i; 
return  c; 

} 


fcomplex  RCmuI(float  x,  fcomplex  a) 
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{ 

fcomplex  c; 
c.r=x*a,r; 
c.i=x*a.i; 
return  c; 

} 

/*  The  error  routine  for  numerical  recipes  routines  */ 
void  nrerror(char  error_text[]) 

{ 

fprintf(stderr,  "run-time  error. .  .\n"); 
fprintf(stderr,  ''%s\n"  ,error_text) ; 
fprintf(stderr,"...now  exiting  to  system... \n"); 
exit(l); 

} 


/*  These  are  utility  functions  for  defining  vectors  and  matrices*/ 
float  *vector(int  nl,int  nh) 

{ 

float  *v; 

v=(float  *)malloc((unsigned)  (nh-nl+l)*sizeof(float)); 
if(!v)  nrerror("allocation  failure  in  vector()"); 
return  v-nl; 

} 

int  *ivector(int  nl,int  nh) 

( 

int  *v; 

v=(int  *)malloc((unsigned)  (nh-nl+1) 

*sizeof(int)); 

if  (!v)  nrerror("allocation  failure  in  ivectorQ'Oi 
return  v-nl; 

} 

void  ffee_vector  (float  *v,int  nl,int  nh) 

{ 
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} 


free((char*)  (v+nl)); 


void  free_ivector  (int  *v,int  nUnt  nh) 

{ 

free((char*)  (v+nl)); 

} 


float  **matrix(int  nrl,int  nrh,int  ncl,int  nch) 

{ 

int  i; 

float  **m; 

m=(float  **)  malloc((unsigned)  (nrh-nrl+l) 
*sizeof(float*)); 

if  (!m)  nrerrorC'allocation  failure  1  in  matrix()"); 
m  -=  nrl; 


for  (i=nrl;i<=nrh;i++)  { 

m[i]=(float*)  nialloc((unsigned) 

(nch-ncl+1)*  sizeof(float)); 
if  .(!m[i])  nrerrorC'allocation  failure  5  in  matrix()"); 
ni[i]  -=  ncl; 

} 

return  m; 

} 

void  free_matrix(float  **m,int  nrl,int  nrh,int  ncl, int  nch) 

{ 

int  i; 

for  (i=nrh;i>=nrl;i++)  free  ((char*) 

(m[i]+ncl)); 
free  ((char*)  (m+nrl)); 

} 
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